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ABSTRACT 

The injection and evolution of turbulence in the intergalactic medium is studied by 
means of mesh-based hydrodynamical simulations, including a subgrid scale (SGS) 
model for small-scale unresolved turbulence. The simulations show that the production 
of turbulence has a different redshift dependence in the intracluster medium (ICM) and 
the warm- hot intergalactic medium (WHIM). We show that turbulence in the ICM is 
produced chiefly by merger-induced shear flows, whereas the production in the WHIM 
is dominated by shock interactions. Secondly, the effect of dynamical pressure support 
on the gravitational contraction has been studied. This turbulent support is stronger 
in the WHIM gas at baryon overdensities 1 ^ 5 < 100, and less relevant for the ICM. 
Although the relative mass fraction of the gas with large vorticity is considerable (52% 
in the ICM), we find that for only about 10% in mass this is dynamically relevant, 
namely not associated to an equally large thermal pressure support. According to 
this result, a significant non-thermal pressure support counteracting the gravitational 
contraction is a localised characteristic in the cosmic flow, rather than a widespread 
feature. 

Key words: Hydrodynamics - Methods: numerical - Turbulence - Cosmology: large- 
scale structure of Universe - Shock wave 



1 INTRODUCTION 

■ Often considered just a by-product of the virialisation mech- 
' anism, the injection of turbulence energy in the intergalac- 
, tic medium (IGM) has been recognised as an interesting 
' process in its own right. Put in a broader context, the in- 
formation about non-thermal phenomena inside and out- 
side galaxy clusters can complement the information derived 
from the hot, X-ray emitting gas. Turbulence, in particular, 
is of utmost importance because of the apparent link to the 
acceleration mechanisms of cosmic rays (CR) and to clus- 
ter diffuse radio emission IIFerrari et al ]|2008l:ICassa"noll2009l : 
iBrunettil l2009l : iBrunetti et all l2009l : ICassano et al] |2010D . 
Turbulent velocity fluctuations are also considered to be 
important for a mplifying magnetic fields i n the intracluster 
medium flCM) llSubramanian et al.ll2006l ) and in the IGM 
||Rvu et al.ll2008l ). 

To date, constraints on the turbulent velocity in clus- 
ters have been d etermined by m easuring resonant scatter- 
ing suppression (jChurazov et al.i i2004 : .Werner et al., ,200 9) 
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and by XMM-Newton spectroscopic observations of clusters 
with a compact core l|Sanders et al.ll2010l . [2OI1I ). Recently, 
the Suzaku satellite made deep X-ray obser vations out to 



the v i rial radius of several clusters pos sible fRcipric h et al 
20091 : iGeorge et al.l |2009|: iBautz et al.l 12009 : .Hoshino et al . 
201ol : lKawaharada et ablboiol ). thanks to its low and stable 



particle background. In most of the published data, the ob- 
served clusters show departures from th e hydrostatic equi- 
libriu m in their outskirts, in one case l|Kawaharada et aU 
the non-thermal support is as large as about 50% of 
the total pressureQ. 

As for the gas outside clusters, the strongest recent ar- 
guments in favour of a significant effect of small-scale turbu- 
lence have been derived from \ow-z O VI observations com- 



^ By 'non-thermal pressure' one indicates generically pressure 
contributions from turbulence, magnetic fields and non-thermal 
particles (cosmic rays); in the following, the latter two contribu- 
tions will not be addressed. Moreover, non-equipartition effects of 
the low-density plasma in the cluster outskirts could mimic a non- 
therm al contribution fsee lWone &: SarazinI [20091 ; iRudd fc Naeail 
l2009i . and references therein). 
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pared with numerical simulations (|Oppenlieimer fc Pavel 
l2009 h. In this work, turbulence had to be added by hand 
as a post-process correction since a subgrid-scale turbulence 
model was unavailable. 

Whereas the theoretical and numerical exploration of 



fair amount of attention (e.g 


, Dolag 


et al. 2005': Vazza et al. 


200f: 


ISubramanian et alJ 


20061: 


Scannaoieco & Briieeen 


2008 


lapichino & Niemever 


l2008l: 


Briiggen & Scannapiecc 


200C 


Vazza et alJ 2009bl: Maier et alJ 20091). the turbulent 



state of the cluster outskirts and of the warm-hot inter- 
galactic medium (WHIM), which is believed to contain a 
significant fraction o f the baryons in the low-z universe 
l|Cen fc Ostrikeij|l999l : IOppenheimer fc Dav^l2008l ). has not 
yet been addressed in a similarly systematic way. The prop- 
erties of turbulence in the cl uster outskir ts have been studied 
recently in simulations by iBurns et all (|2010l ). who found 
that those regions are not in hydrostatic equilibrium, and 
have a substantial turbulent pressur e support. Based on a 
sample of sixteen cluster simulations, iLau et al] (|2009l ) con- 
clude that the support of tu rbulent mot i ons in creases to- 
wards the cluster periphery. IShaw et al.l (|2010| ) find that 
the previous results lead to a significant reduction of the 
Sunyaev-Zel'dovic h power spectrum at angular scales of 
a few arcminutes. IZhu et all l|2010l ) present an instructive 
analysis of the resolved vorticity field and the influence of 
the turbulent pressure in cosmological simulation. Using the 
dyn amical equation for the rate of change of the divergence 
(cf. ISchmidt|[2009bl ). they estimate the effect of turbulent 
pressure on the gravitational contraction of the baryonic gas. 
In this way, the role of gas turbulence in the c lustering pro- 
cess is studied (see also iBonazzola et al.lll99^ ). 

There are several mechanisms that are potentially able 
to stir the baryons and inject turbulence in the fluid. In the 



framework of our simu 
motions in the ICM 


lations we neglect the effects of f 


galaxy 


(iBregman fc DavidI 19891: iKim 


20071: 


iParrish et al.l l2010l: iRuszkowski fc OhI 


20101) and outflows 


from AGN activity ( 


Heinz et al.1 2006 


: Siiacki fc Springell 


l2006l: iBriiEeen et al.l 

shocks thus remain as 


20091). Cluster mergers and curved 
the main stirring agents to be consid- 



ered. 



According to the hierarchical scenario for clustering, a 
halo accretes most of its mass by mergers. In particular, one 
can distinguish between an earlier phase of major mergers, 
where the ratio of the mass of the merging subclumps is 
close to unity, and a subsequent minor merger phase, when 
smaller subhaloes fall into the cluster gravitational well. 

Both mergers phases perturb the cluster medium and 
are thus related t o the injection of turbulence in the 
ICM. According to ISubramanian et al.l (|2006l ). turbulence 
produced in the major merger phase has a large vol- 
ume flUing factor, as expected from even ts which deeply 
stir a nd re a rrange the cluste r stru c ture (iRoettiger et al.l 



19931. 1997 


: Kicker fc SarazinI l200ll 


Paul et al.l 


2011). In case of minor 



of idealised setups has shown, the shear at the bound- 
ary between the ICM and the accreting subcluster trig- 
gers the Kelvin-Helmholtz instability, locally i njecting tur- 
bulence in the wake of the moving subclump s l|Heinz et al.l 
I2OO3I : llapichino et al.ir2008l : lMaier et al.ll2009l ). 

Another stirring mechanism is linked to the baroclinic 
vorticity generations at curved shocks. As known from the- 



ory, vorticity uj is produced where the pressure and density 
gradients are not pa rallel. Taking the curl of the Euler equa- 
tion (tLan dau fc Li fshitz 1959; Kan g et al. | |2007l ). 



du: „ , , Vp X Vp 
_.Vx(.x^)-— ^ 



(1) 



where the second term at the right-hand side is non- 
vanishing at the locations of non-planar shocks. In fllaments 
and cluster outskirts the unprocessed gas is accreted from 
the voids and accelerated towards the growing structures, 
where accretion shocks are formed. The injection of turbu- 
lence is therefore a by-product of the gas accretion at curved 
shocks. 

The main difference in the two distinct mechanisms de- 
scribed above is in the driving: in merger events, turbulence 
is generated by shearing instabilities, whereas at shocks the 
generation is driven by compressional modes. This difference 
is expected to affect the flow features in a quantiflable way 
that will be explored through our numerical simulations. 
Even for a flxed temperature, numerical studies of forced 
supersonic turbulence indicate significant differences in the 
distributions of density and ve l ocity fluctuations dependin; 
on the forcing ("Schmidt et al.' '2008'; 'Peder rath et al.l I2OO; 
ISchmidt et"an [2009: Fodcrrath ot al. 2010). 

It is clear from the above that turbulence is generated 
on essentially all cosmologically relevant scales. Moreover, 
the turbulent WHIM and ICM include a temperature range 
from 10^ — 10® K and corresponding Mach numbers rang- 
ing from 10~^ to transonic values. Adding to the inherent 
difficulties of simulating turbulence even under more sim- 
plified (i.e., homogeneous and isothermal) conditions, these 
complications make the exploration of large-scale structure 
turbulence challenging, to say the least. A convincing sta- 
tistical analysis of the turbulence properties of the IGM on 
numerically resolved scales is infeasible with present codes 
and resources. 

On the other hand, much can be learnt already from 
looking at the magnitude of the production terms and keep- 
ing track of the amount of turbulent kinetic energy on un- 
resolved scales by means of a subgrid-scale (SGS) model. 
The central simplifying assumption here is that turbulence 
can be considered to be statistically isotropic on sufficiently 
small scales. Clearly, this is only a first approximation which 
cannot replace the information gained by increased resolu- 
tion. However, extensive experience with large-eddy simu- 
lations (LES), as this technique is commonly referred to, 
shows that many properties of unresolved turbulence (most 
importantly, those related to transport and dissipation) can 
be captured wi th a certain degree o f confidence (for further 
references, see ISchmidt et al. I I2OO6I : ISchmidtl l2009al ) . Once 
a reliable model for small-scale turbulence in cosmological 
simulations is available, many of the questions listed above 
can already be addressed to some extent. This is the ap- 
proach we will take in t his work. 

In a previous paper (|Maier et al.ll2009l ). we introduced a 
version of LES suitable for adaptive mesh refinement (AMR) 
called FEARLESS. Its main properties will be summarised in 
Section (2] Here, we apply this model to a large-scale struc- 
ture simulation with adiabatic gas dynamics for the first 
time. Considering only the SGS turbulence energy as a probe 
for the production of turbulence for the reasons explained 
above, we find that its evolution indeed differs significantly 
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for the WHIM and ICM phases of the IGM. Although this 
result is limited by the approximate nature of the turbu- 
lence SGS model, we offer an interpretation in terms of dif- 
ferent dominant production mechanisms in these phases and 
present supporting evidence in Section [S] 

In Section |4l we will further elaborate on this approach. 
The turbulence SGS model allows us to predict the contribu- 
tion of the turbulent pressure on the grid scale in addition to 
the effects caused by numerically resolved turbulence. Com- 
paring the WHIM and the ICM, we find that the support 
of the gas against gravitational contraction by turbulence 
is more pronounced at low densities in the WHIM than at 
higher densities in the WHIM and ICM, in which the sup- 
port is mainly thermal. A further important finding which 
will be discussed is that the turbulence-supported gas has a 
fairly low mass and volume fraction. 

We conclude with a summary of our results and sugges- 
tions for future directions in Section [5] 



2 NUMERICAL TOOLS 

This work is based on hydrodynamical simulations of 
the evolution of the cosmic large-scale structure, per- 
formed using the fearless numerical technique (Fluid 
mEch anics with A daptivcly Refined Large Eddy Simula- 
tions; iMajeF^eT^Lji^OOQ.) for simulating intermittent turbu- 
lent flows in clumped media. This tool has been implemented 
on the public release of the grid-based, adaptive mesh reflne- 
ment (AMR) hybrid (N-Body pl us hydrodynamical) code 
ENZO (v. 1.0) (|0'Shea et al.ll2005lfl . 



2.1 Setup of the simulations 

A flat ACDM cosmology is assumed, with Qa = 0.721, S7m = 
0.279, fib = 0.046, h = 0.7, as = 0.817, and n = 0.96. 

The computational box has a side of 100 Mpc and 
is resolved with a root grid of 128'^ cells and 128^ N-body 
particles. The mesh is refined with four additional AMR lev- 
els (refinement factor N — 2), leading to the effective spatial 
resolution of Iaa ~ 48.8 kpc . The force resolution of the 
gravity solver is of the order of 2 x Iaa- The AMR criteria 
are based on bar yon and DM ove rdensity, with overdensity 
factors / = 4 (lapichino fc Nieme vcr 2008). 

The initial redshift of the simulation s is z = 60, and 
the in itial conditions are produced with the lEisenstein fc Hul 
1 19991 ) transfer function. The evolution is then followed to 
2 = 0. Additional physics such as cooling, feedback and 
transport processes is neglected. An ideal equation of state 
was used for the gas, with 7 = 5/3. 



2.2 Subgrid-scale model and the fearless 
approach 

fearless combines AMR with a subgrid scale (SGS) model 
for the unresolved turbulence energy, which encompasses the 



ENZO homepage: |http://lca. ucsd.edu/portal/software/enzo] 



production, the diffusion and the dissipation of kinetic en- 
ergy on subgrid scales (see lSchmidt et al]|2006l ). Details, nu- 
merical tests and applications to the physics of g alaxy clus- 
ters are presented elsewhere (|Maier et al.ll2009l ). and here 
we rec all the main fea t ures of this tool. 

In lSchmidt et all (|2006l ) it is shown how the governing 
equations of a compressible, viscous, self-gravitating fluid 
can be decomposed into a large-scale (resol ved) and a small- 
scale (unresolved) part by exploiting the ICermanol (|l992l ) 
filtering form alism, applied to density- weighted variables 



(Pf) = {p)f^f 



(2) 



(lFavrelll969l ). According to this formalism, once a filtering 
length is set, a variable / can be decomposed in a smoothed 
part {/) and a fiuctuating part /', with {/) varying only on 
scales larger than the filter length scale. A filtered quantity 
/ is thus defined by 

M 
(p) ■ 

In the following, we assume that the filter length scale is 
generically given by the grid scale ^a, i. e., the size of the 
grid cells at any level of refinement. 

By applying this formalism o ne can de rive the filtered 
equations of the fluid dynamics fcf. lSchmidt et al. 2006 ). For 
the sake of conceptual clarity, we do not includ e the cosmo- 
logical expansion here and refer the reader to iMaier et al.l 
(2009) for a complete formulation in co-moving coordinates. 
The resulting equations read 



dt 



{p)Vi + ■7T-Vj{p)'"i 
(JT j 



dt 



dr 



d 



(3) 



(4) 



+ {p)9i - T^-rC 



drj 



V^,Vj) 



dt 



+ 



d 



+ + <:) - ■Ci— f (u,,uj) 



+ {p)'"i9i 



d . 
r(wj,eint) 



(5) 



dr 



where p{ri,t) is the baryon density, Vi{ri,t) are the velocity 
components and e{ri,t) is the total speciflc energy, p the 
pressure, gi the gravitational acceleration and a'ij the vis- 
cous stress tensor. In the above equations, the generalised 
moments of arbitrary quantities / and g are given by 



(6) 



The flltering of energy leads to the definition of a total 
resolved energy eres = eint + 1/2 ViVi,w\iexe eint is the internal 
energy, and the term 1/2 ViVi is the resolved kinetic energy. 
On the other hand, the filtered kinetic energy ekin includes 
also an unresolved contribution, expressed by a second-order 
moment of the velocity field T{vi,Vj), the turbulent stress 
tensor: 



ekin = ^WiWi + '^T{Vi,Vj)/{p) 



(7) 



As in lGermand l|l992l ) we identify the trace of T{vi,Vj)/ {p) 
with the square of the SGS turbulence velocity q, so that we 
define the SGS turbulence energy as 



et = ^q^ ■- ^T{vi,Vi)/{p). 



(8) 
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Since the trace of T{vi,Vj) can be added to the thermal 
pressure (p) in the fihered momentum equation Q, this 
identity immediately implies that 



Pt = 3(p)et 



(9) 



is the turbulent pressure associated with the turbulent veloc- 
ity fluctuations on length scales smaller than the grid scale 
Ia. 

The governing equation of et , as derived by the filtering 
of the equations of fluid dynamics, is 

lr(p)et + 7^".{p)et= D + E + r-(p)(A + e) , (10) 



dt 



where the terms at the right-hand side have to be explicitly 
defined as a function of large-scale filtered quantities and 
of et. Their definitions (the so-called closures) represent the 
turbulence SGS model. The physical interpretation and the 
expressions of the terms in equation pop are the following: 

• T> represents the diffusion of SGS turbulence energy. 
Its expression is based on the gradient-diffusion hypothesis 
l|Sagautll200'^ ) 



d d 
2? = -^C-D{p}lAq'^-^q 
art dri 



(11) 



with Ct> = 0.4, a s inferred by numerical experiments 
l|Schmidt et al.|[2006l ). and I a is the cutoff scale. 

• E is the production term, i.e. the term which accounts 
for the flux of kinetic energy from the resolved to the SGS 
component. The production terms arises in equation (jlOp 
as the contraction of the turbulent stress tensor and the 
Jacobian of the resolved velocity: 



E = 



(12) 



The turbulence stresses f (ui, Vj) are given by the commonly 
used eddy- visc osity clo sure. In the comp ressible formulation 
ijSchmidt eral., 2006: Maier et al.ll2009f l. this closure reads 



-{v^,Vj) = -2{p)CjAqS*j + -Sij{p)q^ 



(13) 



where 5**, is the trace-free part of the rate-of-strain tensor, 



_1 f dv, dvj 



(14) 



• r accounts for small-scale gravitational effects, and is 
neglected in our implementation. As further simplifications 
in the model, the infiuence of the viscous stress tensor {o-'ij) 
in equations ^ and ((Sj is neglected (which is well justified 
for high Reynolds numbers), and so is also the SGS transport 
of internal energy, given by the divergence of f(wj,eint) in 
equation ([5}; 

• e is the term for SGS turbulence dissipation into inter- 
nal energy, treated as an effect acting only at subgrid scale. 
This term is therefore added to the numerical dissipation of 
the code, but is m odelled in a physically motivated way (see 
iMaier et aLll2009l for i t s role in the physics of the ICM). Its 
closure follows ISarkaij l| 19921 ): 



af-(i + QiM2) 

(A 



(15) 



where Q = 0.5 (|Sagautll200l l, ai = 0.5 and Mt = q/c^ is 
the turbulent Mach number, where Cs is the speed of sound. 



• A is the pressure dilatation term, which models the ef- 
fect of unresolved pressure fluctuations and, effectively, is 
an exchange term between et an d eint. Like fo r the previ- 
ous term, a closur e suggested by ISarkail l| 19921 ) is adopted 
(jMaier et al.ll2009l) . 

An enlightening graphical representation of the SGS 
terms and the i r role in the energy budget is given by 
llapichino et"al] (|2010l . cf. iMaieiiliooj ). 

The innovation of the fearless approach is a consistent 
treatment of the interchange between the kinetic energy on 
resolved length scales and the SGS turbulence energy for 
a cutoff length that varies in space and time. Apart from 
the energy flux through the turbulent cascade, this involves 
the increase of resolved kinetic energy at the cost of SGS 
turbulence energy if refined subgrids are inserted and, thus, 
the numerical cutoff length scale decreases. 

The i nterplay between the AM R of enzo (based on the 
method of iBerger fc Colellal IT989I ) and the turbulence SGS 
model in fearless exploits an additional assumption of Kol- 
mogorov sca ling of the turbulent energy l|Kolmogorovlll94l] : 
lFrisch|[l995l ). Considering two AMR levels with spatial res- 
olutions Za,; and Iaj, it means that, at grid refinement (or 
derefinement) , the SGS turbulence energies are statistically 
related by 



et,» 



^A.i 



2/3 



(16) 



When a region is refined, a new "grid" is created, i.e. the 
refined grid patch is handled as a single AMR object. In this 
new fine grid the values of the hydrodynamical variables are 
interpolated from the coarse grid, but the SGS turbulence 
energy is scaled according to equation (|16p . and the veloci- 
ties are corrected such that the sum of resolved energy and 
turbulent energy remains conserved. An opposi te procedure 
applie s to grid derefinement, as described by iMaier et al.l 
(l2009h . 

An important preliminary check concerns the range of 
applicability of the turbulence SGS model, which is devised 
to the study of, at most, moderately compressible flows. 
This limitation is incorporated in the code as a safeguarding 
mechanism for the value of the turbulent Mach number. In 
order to prevent numerical ins tabilities, a thresho ld is set at 
Aft, max ~ V^, as motivated in lMaier et al.l l|2009h . 

The issue of the applicability of the turbulence SGS 
model for different baryon temperatures is summarised in 
Fig. [T] Most of the cold gas (T < lO" K) lies in the plot at 
Alt ~ \/2: for this baryon phase the low sound speed makes 
the gas motions very supersonic. The turbulence production 
is therefore unphysical for this gas, because the implemented 
SGS model is not suitable for its study, and also because the 
cold and rarefied medium is poorly resolved by our refine- 
ment criteria. Anyway, the study of this cold gas phase is 
physically not well posed in simulations without UV back- 
ground heating. 

In order to circumvent this shortcoming, we will limit 
our analysis to the gas with T > 10^ K. The threshold on the 
turbulent Mach number does not affect the gas above this 
temperature, where Mt is typically below unity. The SGS 
turbulent energy is thus not a dominant c omponent in t he 
energy budget, as already discussed in Maier et al.l l|2009l ). 



© 2011 RAS, MNRAS OOO.mflil 



Turbulence production and pressure support in the IGM 5 



1.6 




10^ lO'^ 10^ lo'" 10^ lO*" lo' 10^ 10^ 



TemjDeratiire |K} 

Figure 1. Two-dimensional mass distribution function of tiie tur- 
bulent Mach number Mt as a function of the gas temperature, at 
the redshift z = 0. The mass is coded according to the colour bar 
on the right. The sharp horizontal line at Mt = 1.41 is caused 
by the safeguarding mechanism in FEARLESS. A second horizontal 
hue, visible at Mt = 1.41 X 2-1/^ = 1.12 and T < 10* K, is an 
artifact produced by the interaction of the threshold on Mt and 
the refinement at the first AMR level. 

3 EVOLUTION OF THE WHIM AND THE ICM 

3.1 Baryon phases 

The main results of this work make use of a distinction of 
the baryons in two phases. This distinction is somewhat ar- 
bitrary, and indeed the mass distribution function in Fig. [2] 
shows that the cosmic gas is characterised by a continuum of 
gas states in temperature and density. The problem is gener- 
ally addressed in the li terature by using a thr eshold based on 
gas temperature (e. g ..lCen fc Ost rikerl^ 19991') and/or density 
IVazza et all l2009al : ISkillman et al. 2008). Ideally, a phys- 
ically motivated distinction should be done on a dynami- 
cal basis, between gas belonging to virialised structures or 
not. This criterion would be computationally too demand- 
ing, and is not used in our analysis. 

In the following we will make a distinction based on the 
baryon overdensity of the gas 5 = p/(f2bPcr), where pcr = 
3Hq/{8ttG) (l + z)^ is the critical density at redshift z. The 
gas will be labelled as 'WHIM' if 5 < 10^ and as 'ICM' if 5 > 
10^. As discussed at the end of Section [2.21 for both phases 
the additional constraint T > 10* K is imposed. The former 
baryon phase is mostly to be found in filaments and in outer 
halo atmospheres, whereas the latter is in the potential wells 
of groups and clusters (cf. Fig. [3]). From the definition of the 
AMR criteria in our simulation, it follows that the ICM gas 
is always resolved at the highest AMR level I a ~ 4, whereas 
the mass-weighted average refinement level for the WHIM 
is 2.6. 

3.2 Turbulent energy 

A first overview of the simulation data a.t z — (Fig. ^ 
already shows a match of the locations where both the in- 
ternal and the SGS turbulent energies are large. This indi- 




10'' 10"' 10° 10^ 10^ 10^ 10" 10^ 

Baiyon Overdensity 

Figure 2. Two-dimensional mass distribution function of the gas 
temperature as a function of the baryon overdensity, at z = 0. The 
mass is coded according to the colour bar on the right. 
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Figure 4. Time evolution of the mass-weighted averages of spe- 
cific internal (ei^t, solid lines) and SGS turbulent (et, dotted lines) 
energies, for the two baryon phases under investigation. The two 
lines in the upper part of the plot refer to the ICM, and the other 
two to the WHIM. The lines are scaled according to the factors 
in the legends, in order to be accommodated in the same plot. 

cation is consistent with the close link occurring between 
gravitational collapse, virialisation and injection of turbu- 
lent energy during cosmological structure formation. In this 
Section, we use the SGS turbulent energy (measuring tur- 
bulent velocity fluctuation at the numerical cutoff scale) in 
relation to the internal energy (measuring the temperature) 
of the gas as a diagnostic of the properties of turbulence. Al- 
though turbulence is produced on length scales larger than 
the grid cutoff scale, the local energy injection is imprinted 
on the SGS turbulent energy because of the energy transport 
through the turbulent cascade. 

In Fig. [4] the temporal evolution of the mass- weighted 
average of the internal (eint) and SGS turbulent (et) speciflc 
energies, both for the WHIM and ICM phase, is reported 
in more detail. From et one can derive the average SGS 
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Figure 3. The panels show projections of a cube with the side of 20 Mpc h~^, extracted from the computational domain of the FEARLESS 
run, at 2; = 0. Baryon overdensity is shown at the left-hand side, internal energy in the central panel, and the SGS turbulent energy at 
the right-hand side, with density contours overlayed, and the corresponding colour bar under the panels. 



turbulent velocities for the WHIM and ICM: they are 59 
and 76 km s~^, respectively, at z = 0. The corresponding 
average turbulent Mach numbers are 0.18 and 0.14 . 

The key feature to be noticed in Fig. |4] is the different 
trend for the evolution of et with time: for the WHIM phase, 
et increases steadily to 2: = 0, whilst in the ICM it reaches 
a peak between z ~ 1.0 and 0.65, and then decreases. In 
the following we speculate that this different evolution is re- 
lated to the mechanisms of turbulence generation in the two 
baryon phases. More specifically, the time evolution of et in 
ICM and WHIM can be interpreted as turbulence produc- 
tion by mergers in the ICM, and by shock interactions in 
the WHIM. 

There is supporting evidence corroborating our hypoth- 
esis. The gas in the ICM belongs to collapsed structures, 
which experienced merger episodes during their evolution. 
These merger events are related with the stirring of the 
baryons and the subsequent injec tion of turbulence (e.g., 
iMaier et ai]|2009l : [Paul et aLlbOllT ). It is therefore not sur- 
prising that the maximum of et in the ICM is consistent 
with the formation time (defined by the major merger phase) 
in the evolution of clusters in the mass range 10^^ Mq < 
M < 10^* MfT), as inferred by analy t ical models of hierarchi- 
cal clustering jLacev fc CopllQQl ISheth fc TormenI |2004 
and Fig. 4 of lCiocoh et al.ll2007l ). At later times, the level 
of SGS turbulence decreases, although the local injection 
in minor mergers likely cont ributes to slow down its decay 
(cf ISubramanian et al.ll200d ). 



|2000| . or Fig. 2 of lSkillman et al.ll200d^T ). This hints towards 
a link between the SGS turbulence energy production in 
the WHIM and the gas accre tion on shocks associa ted with 
growing structures. Recently. ICavaliere et al.1 (|201ll ) investi- 
gated thoroughly the injection of turbulence in cluster out- 
skirts (at the accretion shocks) with analytical calculations, 
finding that the turbulent support increases from z ~ 0.5. 
In their analysis, this is due to the shock weakening, the 
decrease of accretion rates on clusters and the decrease of 
the gas infall speed at low redshift. The evolution of et for 
the WHIM at redshift aroun d 0.5 in Fig. [3] is in qu alitative 
agreement with the model of ICavaliere et al] (|201lh . 

In order to verify the resolution insensitivity of our re- 
sults, we analysed the inn ermost part of the fearless clus- 
ter simulation discussed in lMaier et al.l (|2009l '). In that setup, 
a computational box with a side of 128 Mpc is simulated 
with a root grid resolution of 128^ cells and 128"^ N-body 
particles, but in a small cube with a side of 32 Mpc , an 
additional static grid is nested and seven AMR levels are al- 
lowed, such that the local root grid resolution is equivalent 
to 256'^ for both the mesh and the N-body particles, and 
the effective spatial resolution is 7.8 kpc . This small 
volume is centred on a growing cluster and therefore is not 
representative of a random realisation of the cosmological 
initial conditions, but interestingly in this region the time 
evolution of et and eint is equal to that shown in Fig. [4] Our 
results therefore look robust with respect to an increase of 
spatial and force resolution. 



As for the WHIM, the injection of turbulent energy 
at curved shocks is obviously related to the features of 
the gas accretion on filaments and haloes i.e., more specifi- 
cally, on the amount of kinetic energy processed by the ex- 
ternal shocks. In fact, a qualitative similarity can be no- 
ticed between the evolution of et and the flux of kinetic 
energy through external shocks (Fig. 10 of iMiniati et al.l 



^ The kinetic energy flux in|Miniati ct al.' ('2000') is averaged over 
all Mach numbers, whereas in Skillman et al. (2008) it is sorted 
in Mach number, with results grossly comparable with the former 
study. 
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Figure 5. Evolution of the mass- weighted compressive ratio res 
for the gas in the ICM (solid line) and WHIM (dotted line) tem- 
perature phase. 
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Figure 6. Two-dimensional distribution function of the average 
compressive ratio res in the T — 5 plane, at z = 0. 

3.3 Compressive ratio 

Our investigation, started from the analysis of the subgrid 
turbulent energy in Fig. [l] is furthermore supported by the 
structure of the velocity field at resolved sca les. To this aim 
we d e fine the small-scale c ompressive ratio (jKida fc Orszad 
ll990l : ISchmidt et al.|[2009l) : 

r (17) 

where {cF) and (o;^) are the averages of the squares of 
the divergence and the vorticity of the velocity field, re- 
spectively. This ratio quantifies the relative importance of 
compressional and so lenoidal modes in a fiow. As shown by 
ISchmidt et al] (|2009l ). the values of Tcs tend to be higher for 
compressively-driven turbulence. 

The evolution of res for the two baryon phases is re- 
ported in Fig. [5l As expected from the previous consider- 
ations, the gas in the WHIM phase has a larger value of 




Figure 7. Same projection volume as the panels in Fig. [S] but 
the compressive ratio res is shown. Baryon density contours are 
overlayed. 



the compressive ratio throughout the simulation, with re- 
spect to the ICM phase, indicating a higher contribution 
from compressional modes in filaments and cluster periph- 
eries. A more detailed analysis on the T — 5 plane at 2: = 
(Fig. |6]) shows that the compressive ratio is low at high den- 
sities and temperatures (the ICM), while it is significantly 
higher elsewhere. The compressive ratio is particularly high 
also at the extrema of the overdensity distribution, resulting 
either from strong rarefactions or compressions. 

A similar conclusion can be drawn from the visual in- 
spection of the projection in Fig. [7] res is generally lower 
in clusters, except for localised regions (for example, in the 
cluster at the centre of the projected volume), likely to be 
associated with weak shocks in the ICM. 

3.4 Energy comparison with the adiabatic run 

The role of the turbulence SGS model in our large-scale 
structure simulations and the consistency of the energy bud- 
get in this framework is further investigated with a compari- 
son between the simulation using fearless and an adiabatic 
reference run. In these simulations the flow is mostly sub- 
sonic, and the energy content of the SGS turbulence is glob- 
ally almost negligible (|Maier et al.ll2009l ). Quantitatively, it 
means that, for every baryon phase, the sum of kinetic and 
thermal energy should be approximately equal in the adia- 
batic and the fearless run (in the latter case, the sum is 
extended to the SGS turbulent energy). In fact, good agree- 
ment (mostly within 2%) is found on the global gas proper- 
ties of WHIM and ICM (mass fractions, energy content and 
their time evolution). 

As a more sensitive diagnostic for the detailed energy 
budget in the ICM and WHIM, for the different phases we 
study the evolution of the total energies 

E^^ ek,«p.Vi, (18) 

ijphase 

where ek,i is the value of the specific energy ek in the cell i 
(the index 'k' refers to eint or et; cf. Fig. |4]) and piVi is the 
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Figure 8. Temporal evolution of the quantities Et (solid lines) 
and ASint (dotted lines), defined in the text, for ICM and WHIM. 



baryon mass in the cell i. The sum is performed on the cells 
belonging to a same baryon phase, either WHIM or ICM. 

In Fig. [8] we compare Et with AEint, the difference of 
Eint between the adiabatic and fearless simulations. In 
both phases, Eint (not in the plot) is 40 times or more larger 
than Et. 

We observe that Et and Ai?int are of the same order 
of magnitude during the simulation, both for the ICM and 
WHIM, indicating that the SGS turbulent energy acts as an 
energy buffer between the resolved and unresolved scales. In 
other words, the global decrease of Eint in the fearless run 
is partly balanced by Et, so that the global energy budget 
is nearly unaffected. 

Further physical interpretations of this energy budget 
are difficult, because of the high complexity o f the flow in 
cosmological simulations; we refer the reader to'M aier et al.l 
(,200R ) to more tests in simplified setups. Another caveat, 
however, is that the turbulence SGS model is pushed to 
its limit in the WHIM because of the relatively high Mach 
numbers in the fiow. For this reason, the result in the WHIM 
has to be confirmed with a SGS model that does not suffer 
from such constraints (see lSchmidt fc Federrathll201ll ). 



4 THERMAL AND TURBULENT PRESSURE 
SUPPORT 

Using data from a cosmological hydrodynamic simulation, 
IZhu et all (|2010l ) present an in-depth analysis of the vortic- 
ity and divergence fields in the intergalactic medium. The 
rationale behind their analysis is similar to ours, except that 
they infer turbulence properties from the derivative of the 
resolved velocity field. Moreover, they consider dynamical 
equations for the modulus of the vorticity and the diver- 
gence. Of particular importance is the rate of change of the 
divergence, which is generalised to a co-moving coordinate 
system: 



Dt a 



2 P P 



1 



A-KGip + Pdm) 2~^" 



(19) 



where a is the time-dependent cosmological scale factor, 
D/Df = d/dt + a~^v ■ 'V is the material derivative in co- 
moving coordinates, G is the gravitational constant, pdm is 
the local dark matter density, Om is the cosmological mean 
density parameter of baryonic and dark matter, respectively, 
a nd H = a/a is th e Hubble parameter. This is the same as 
in lZhu et alj (|20ld ) (their equation 3), with only the gravity 
terms (in the second line of equation [19} slightly rearranged, 
and noting that, for the components of the rate of strain ten- 
sor, 2SijSij = l^p. 

The advantage of the filtering approach outlined in Sec- 
tion [22] is that we can easily include SGS terms, in particu- 
lar, turbulent pressure terms. The SGS model described in 
Section [2.21 allows for a direct computation of the turbulent 
pressure that is associated with the grid scale: pt = 2/3 pet. 
Since the divergence equation is derived from the momentum 
equation, in which the turbulent pressure is simply added to 
the thermal pressure, it follows that the filtered version of 
the divergence equation is readily obtained from equation 
([19} by substituting p with p + pt everywhere: 



Di a 



AnG{p + pdn 



-Pt, 



(20) 



and by considering filtered quantities (we dropped the hats 
for brevity). The trace-free part of the turbulence stress ten- 
sor is neglected in the above equation. The expression on the 
right-hand specifies the net negative compression rate of a 
fiuid parcel. The self-gravity term on the very right stems 
from the Poisson equation for the gravitational potential and 
tends to decrease the divergence. 

To understand the meaning of the various terms in 
equation (|20p . we consider different limiting cases: 

(i) Incompressible limit: The fiuctuations of the density 
with respect to the mean density vanish and d — 0. Thus, 



1 , 2 

2 ^ 



-\Sf) = -V\p + pt) 
P 



(21) 



(ii) Infinite resolution {I a — )■ 0): In this case, pt vanishes, 
and equation (|19p is obtained. 

(iii) Global filtering (Za ~ L, where L is the integral 
l ength scale for turbulence injection): The filter formalism 
(jGermandlfgg^ ) encompasses the limit of a statistical the- 
ory. If flow structure is smoothed over the largest scales, i. e., 
the size of galaxy clusters, the velocity derivative becomes 
negligible, and the effect of turbulence is entirely given by 
the turbulent pressure: 



-iv^(p + Pt) + ^Vp ■ V{p + pt) 
p p^ 



gj^2 

4:TtG{p + Pdm) 



(22) 



Hd, 



Neglecting the effects of pressure gradients that are un- 
aligned with the density gradients and comparing the limit- 
ing cases (ii) and (iii), we see that the term 1/2 p{uj^ — \S\^) in 
a fully resolved simulation is equivalent to — V'^pt if the fiow 
is filtered on the largest scale of the system. In a large eddy 
simulation, we have an intermediate case, where part of the 
effect of turbulence is captured by the vorticity and the rate 
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of strain of the resolved flow, while the turbulent pressure 
at the grid scale accounts for numerically unresolved turbu- 
lence. If w > \S\, numerically resolved turbulence counter- 
acts the gravitational contraction of the gas. The turbulent 
pressure of unresolved velocity fluctuations counteracts self- 
gravity if V^pt < 0, respectively. The relative contribution 
of pt depends on the grid scale. 

It is important that, by its very definition, the 

turbulent pressure is a scale-depe n dent quantity 

l|Schmidt fc FederrathI HoTl). IZhu et al.l (|2010l ) investi- 
gate the scale-dependence of the turbulent pressure by 
integrating the spectrum of the kinetic energy density for 
all wave numbers greater than then a certain wave number 
(corresponding to a particular length scale). Since the 
resulting turbulent pressure spectrum is rather flat, no 
clear distinction is made between the integral turbulent 
pressure of the resolved flow and the turbulent pressure of 
velocity fluctuations below the grid scale. The advantage 
of our approach is that we can investigate both resolved 
turbulence and SGS turbulence effects. In Fig. (9] we show 
the mass- weighted correlation diagrams of 1/2 p(tJ^ — \S\^) 
and -VV- Both for the WHIM and the ICM, these 
quantities are roughly correlated. This is expected, because 
SGS turbulence is produced by the interactions with 
turbulent velocity fluctuations on the smallest resolved 
length scales, which are probed by uj and However, the 
non-local nature of the SGS turbulence energy (see equation 
nop , implies that there is no simple algebraic relationship 
between the resolved and unresolved turbulent pressures. 
This becomes manifest in the large scatter of the correlation 
diagrams. Consequently, the SGS model is essential for 
the computation of the support by the turbulent pressure, 
—\7^Pt- Compared to a given value of the resolved turbulent 
pressure (corresponding to a horizontal cut through the 
two-dimensional distribution), — V^pt is typically an order 
of magnitude smaller. Locally, however, the contribution 
from SGS turbulence can become comparable to or even 
exceed the resolved contribution, with the exception of 
extremely strong turbulence intensity in the upper right 
corner of the distribution (1/2 p{uj^ - \Sf) > 10"^" and 
10"'^^, in arbitrary units, for panels (a) and (b) in Fig. [51 
respectively) . 

To analyse the effects of intense vorticity rela tive to the 
suppo rt of the gas due to its the thermal pressure. IZhu et al.l 
l|2010ll calculate the ratio of 1/2 p{uj^ - \Sf) to -V^pQ Ac- 
cording to the above discussion, the full effect of turbulence 
is given by the sum of the vorticity term and the negative 
Laplacian of the turbulent pressure on the grid scale. To esti- 
mate the importance of the turbulent relative to the thermal 
pressure support, we plot in the panels of Fig. [TD]the average 
in the T — S plane and the two-dimensional mass distribution 
function of the ratio 

ip(.;^-|S|^)-VV _ (23) 



Table 1. Mass (second column) and volume fractions (third col- 
umn) of WHIM and ICM gas, selected according to the definitions 
(first column) introduced in Section U] at z = 0. 



rtp 

subject to the constraints 
1 



V2p 



p{iu^ -\Sf)-V'^Pt>0 and -V^p>0 (24) 



•* However , V^p appears with a positive sign in figure 10 of 
IZhu et al ] l|2010l l. It is not clear whether this is an error in the 
labelling or an inconsistency of their calculation. 





mass fraction 


volume fraction 


f^WHIM 


0.315 


1.83 X 10- 


-2 


-PL;, WHIM 


0.287 


0.221 




^rtp.WHIM 


7.74 X 10^2 


5.95 X 10- 


-2 


f^ICM 


7.98 X 10-2 


3.53 X 10- 


-5 




0.523 


0.496 




■^rtp,ICM 


0.128 


0.117 





The distribution in Fig. llOf a) suggests higher values of 
rtp for overdensities S ~ 10. To be able to discern this trend 
more clearly, in Fig. llOf b) the mass distribution function of 
rtp as a function of 5 is shown. A bimodal mass distribution 
becomes apparent, where the low-density peak is a feature of 
the WHIM, and the high-density peak is related to the ICM. 
Besides fluctuations, the average of rtp decreases for increas- 
ing baryon overdensity. Clearly, the turbulent support (re- 
solved plus unresolved) is important for the low-density gas 
in the WHIM (1 < 5 < 100). For this gas, the average rtp is 
larger than for the WHIM at higher overdensity, and for the 
ICM. This fact can be linked with the energy evolution seen 
in Fig. 21 the driving of large-scale turbulence in the ICM at 
2 = has passed its maximum and is declining, while turbu- 
lence production in the WHIM is just saturating. Moreover, 
at low redshift et is a larger fraction of eint in the WHIM, 
rather than in the ICM (according to Mt for the two phases, 
cf. Section [3]2]). In the ICM we observe a declining trend of 
rtp towards large overdensities, meaning that the support 
for the densest gas is, on the aver age, mainly ther mal. 

A more precise comparison to lZhu et al.l (|2010l . their fig- 
ure 10) is not straightforward, because we use mass distribu- 
tion functions in place of their scatter plot, and we include 
the contribution of SGS turbulence to the total pressure. 
Moreover, since the WHIM and the ICM are not clearly sep- 
arated in their plot, the trend with the density is obscured. 
Anyway, the a verage values in F ig. [TOl are qualitatively sim- 
ilar to those of lZhu et al.l l|2010l l. making us confident of the 
consistency of these two analyses. 

It is very instructive to consider the relations in equa- 
tion (|24p separately. The first one defines the locations where 
the flow is turbulence-supported, regardless of the thermal 
pressure support on the divergence equation. We notice that 
this relation selects regions where the vorticity of the flow 
uj is large, thus recalling one of the outstanding features of 
turbulence. 

In Table [1] we report the mass fractions of the WHIM 
and ICM gas JIwhim and flicM , and the fractions of the com- 
putational volume occupied by this gas. Then, with Fui^wmm 
we define the fraction of fiwHiM where the first relation in 
equation (|24p is fulfilled, and the same for the ICM gas. 
Both for the WHIM and ICM, is a considerable frac- 
tion of the gas, in mass and volume. Moreover, one can 
see that the ratio between the volume and mass fraction 
is much smaller for Q than for _F„, because the 'cj-selected' 
gas is much less clumped than the whole density-selected 
gas. This is in agreement with Fig. 1101 where we observed 
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Figure 9. Mass-weighted correlation diagrams of the Laplacians of the "resolved turbulent pressure", 1/2 - |5|2)^ and the turbulent 
pressure on subgrid scales, — V^pt, for the WHIM (a) and the ICM (b) at redshift 2: = 0. In both panels, the diagonal line marks the 
location of the equality 1/2 p{u)'^ — \S\'^) = — V^pt- 
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Barvon Overdonsity Baiyon Ovcrdoiisity 



Figure 10. Panel (a): two-dimensional distribution function of the average pressure ratio rtp, defined in equations I I23I I and 1124 l l. in 
the T — 5 plane, , at 2: = 0. Panel (b): two-dimensional mass distribution function of rtp as a function of baryon ovcrdensity. The 
mass- weighted average is overlayed. 



that the turbulent support is smaller for the densest gas, in 
both phases. 

We also observe that the volume fraction of 
Fui.iCM is rather similar to the volume filling factor 
for turbulence in the I CM estimated analytically by 
ISubramanian et al. ||2006|). an d com puted in the simula- 
tions of llapichino fc Niemeverl l|2008l V The latter analysis 
was based on a simpUfied definition for characterising the 
turbulent flow based on its vorticity, therefore it is not sur- 
prising that it compares well with the present study based 
on the turbulent pressure support. 

As interesting as it may be, the information brought by 
Fu, has to be complemented: given the tight link between 



turbulence injection and thermalisation, it is meaningful to 
consider the role of the turbulence only in combination with 
the thermal component of the energy budget. To this aim, in 
equation (|23|) we defined the ratio rtp between the dynami- 
cal and the thermal terms opposing gravitational collapse in 
the divergence equation pUll . Conservatively, we will assume 
that for rtp > 0.1 the role of the dynamical support in equa- 
tion (|20|) becomes non-negligible. According to this thresh- 
old, let frtp.wHiM be defined as the fraction of fiwHiM where 
both relations of equation p4p are fulfilled, and rtp > 0.1. 
An analogous quantity is defined for the ICM, and both are 
reported in Table [T] in mass and volume fractions at ^ = 0. 
Although the mass and volume fractions expressed 



© 2011 RAS, MNRAS 000, [THM] 



Turbulence production and pressure support in the IGM 11 



by Fui are relatively large, the turbulent support is non- 
negligible only for a smaller subset of these cells. Frjp for 
the ICM is of the order of 10% of the global mass and vol- 
ume, and somewhat smaller for the WHIM. This result sug- 
gests that a significant non-thermal pressure support, coun- 
teracting gravitational contraction, is a local behaviour of 
the cosmic flow, rather than a widespread feature. 



5 CONCLUSIONS 

In this work we studied the evolution of the energy budget 
of the ICM and WHIM in mesh-based hydrodynamical cos- 
mological simulations. Besides the internal and the resolved 
kinetic energy, the fearless method combining AMR and 
LES allows us to also follow the evolution of the unresolved, 
SGS turb ulent energy, define d according to the model de- 
scribed in iMaier et al.l l|2009l ). Since an energy cascade sets 
in from the integral length scale down to the (unresolved) 
Kolmogorov scale, the SGS turbulence contains informa- 
tion on turbulent injection and evolution at larger, resolved 
scales, with the advantage of being computed cell-wise and 
thus easily accessible without further post-processing the re- 
solved gas velocity data. 

A first result of this work is the indication of a pro- 
duction of turbulence with different properties in different 
baryon phases (Fig. 3]). In the ICM, the SGS turbulent en- 
ergy peaks approximately at the expected formation redshift 
of the haloes with mass 10^'^ - lO'^* Mq, indicating a tur- 
bulence driving mechanism associated with merger events. 
Indeed, the compressive ratio of the ICM baryon phase is 
relatively low (Fig.[S|, as expected in a flow where the driv- 
ing mechanism is dominated by shearing motions. 

In the WHIM phase, the SGS turbulent energy grows 
more smoothly with time, hinting towards a different pro- 
duction mechanism. It is straightforward to call for the 
role of shocks (in particular, external shocks) in this pro- 
cess, because they enclose the WHIM gas in filaments and 
outer cluster regions. Interestingly, the flux of kinetic energy 
throug h t he external shocks, as simulated by iMiniati et al.l 
l|200(tl and lSkillman et all (|2008l ). closely resembles the tem- 
poral evolution of et . A similar trend has bee n also predicted 
analyt ically in the cluster peripheries by ^ Cavaliere et al.l 
l|201ll l. Further evidence is provided by the compressive ra- 
tio for the WHIM, which is larger than the values found for 
the ICM and thus indicates a flow dominated by a compres- 
sive driving mechanism. The energy content of SGS turbu- 
lence is larger for the ICM phase (Fig. 3}, but the relative 
importance with respect to the internal energy is larger for 
the WHIM, as testified also from the slightly larger turbu- 
lent Mach number (Section 13. 2p . 

Some cautionary remarks are needed for a correct in- 
terpretation of the results. First of all, this bimodality 
is not unambiguous. As known, weak shocks are ubiq- 
uitous in the ICM fMiniati et al] I2OO0I: IRvu et all l2003l: 



Pfrornmer et al.l I2OO6I : ISkillman et all 1200^ : IVazza et al.l 
2009al ). thus adding a compressive contribution in that 



baryon phase, and similarly sm all clumps move also along 
the filaments (|Dolag et al.ll200^ . contributing to solenoidal 
driving modes in the WHIM. Nevertheless, the mixing of 
compressive and solenoidal modes is dominated by the for- 
mer for the WHIM and by the latter for the ICM. It would 



be interesting to explore if and how further stirring mecha- 
nis ms like AGN outfl ows change the simulated scenario. 

IZhu et al.1 (|2010l ') investigate the effects of the turbu- 
lent pressure on the rate of change of the divergence. This 
quantity turned out to be important for the clustering of the 
baryonic gas in the gravi t ationa l wel ls of the dark matter . 
FoUowing llapichino et al] (|2008l ) and lSchmidt et al.1 l|2009D . 
we propose to utilise the negative growth rate of the diver- 
gence as control variable for adaptive mesh refinement in 
cosmological simulations. 

The study of the divergence equation (|20[) has been ex- 
tended with the SGS turbulence modelling, finding that: 

• In general, the contribution of unresolved pressure to 
equation H20|) is not a dominant term, although locally it 
can become relevant (Fig. 

• The turbulent support (defined by the first relation in 
equation I24|) is largest for the WHIM gas in the overdensity 
range 1 < 5 < 100, and tends to be small for the ICM gas 

(5>10^); 

• A non-negligible mass fraction of the WHIM (28.7%) 
and ICM (52.3%) is characterised by a large vorticity, ac- 
cording to the criterion given by the first relation in equation 
(|24l) . However, in most of this gas the thermal support is not 
smaller than the dynamical one, so that the mass fractions 
of the gas where rtp > 0.1 are only 7.7% and 12.8% of the 
WHIM and ICM at 2: = 0, respectively. 

For completeness, we note that at 2: = in our simula- 
tion about 31% of the gas is in the WHIM phase, and 8% in 
the ICM phase: most of the turbulence-supported gas is thus 
in the WHIM phase, although the ICM gas is more tightly 
related to observables. The emerging picture is that of tur- 
bulence as a low-redshift (mostly z < 0.5) feature in the 
WHIM, while in the ICM the main stirring epoch is slightly 
earlier. The total mass fraction of baryons with T > 10^ K 
where the first relation in equation (|24|l holds is 13.2%, in 
good a greement with the previous investigation of lZhu et al.l 
(|2010t ). 

The fact that Fu, > Fr^^ means that the turbulence- 
supported gas is a substantial fraction of the cosmic baryons, 
but the flow producing significant non-thermal support fill 
only small fractions of space. Although this result reflects a 
typical feature of turbulence and can be understood by con- 
sidering that turbulence injection is a by-product of ther- 
malisation during the structure growth, this is apparently 
at odds with recent observational and theoretical claims 
of non-thermal p ressure support, especially in the clus- 
ter outskirts (e.g. iKawaharada et al.ll2010l : iLau et aL I I2OO9I : 
ICavaliere et al.ll201lh ." From the presented results, it seems 
that turbulence has a dynamical role only in a volume frac- 
tion of about 6% of the gas with T > 10"" K. One pos- 
sibility is that the spatial resolution of our simulation is 
not adequate to model the flow close to the accretion shock 
of growing clusters, which plays an impor tant role in in- 
jectin g turbulence in the cluster outskirts (jCavaliere et al.l 
l201lh . However, we tested the results in Table [T] in the well- 
resolved, innermost region o f the computational box pre- 
sented in iMaier et al.l l|2009l ). obtaining basically the same 
results (see also Section [3. 2p . 

We observe that the diagnostics used here for the study 
of the turbulent support (pressure Laplacians) are differ- 
ent from other studies, based on ratios between turbulent 
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an d internal enersie s (or, equivalently, pressure ratios). Only 
in iLau et al] ('2009') , as far as we are aware, the analysis 
is complemented by the computation of pressure gradients, 
the very quantitie s that appear in the cluster mass estimate 
l|Rasia et al.ll2006t ). The discrepancy of the definitions could 
explain why turbulence 'support' (in the sense of the pres- 
sure ratio) is found large and widespread in many galaxy 
cluster simulations, although the dynamical role elucidated 
in the present work is much less significant. One could guess 
that diagnostics like the energy or pressure ratios track the 
turbulent gas in a way similar to F^i, rather than JVtp. 

This work should be considered an exploratory study, 
to be complemented by future simulations, possibly includ- 
ing a more sophisticated treatment of additional physics 
(ionisation background, radiative cooling, galactic winds). 
However, the physical processes discussed in this work are 
governed primarily by gravity, thus the role of additional 
ph ysics is not expec ted to change the scenario drastically 
(cf. iKang et al There is also room for improvements 

in the SGS modelling of turbulence, focusing in particular 
on compression- dominated and inhomogeneous turbulence 
l|Schmidtll20"09al ). As mentioned above, a useful complement 
to this study will be a similar analysis of turbulent pressure 
and turbulence support in simulations of single clusters. This 
is left for future work. 

The study of turbulence in many branches of astro- 
physics has been making progress in recent years, both from 
the theoretical and the observational viewpoint. The more 
one goes into the details of this field, the clearer it becomes 
that a simplistic way of approaching turbulen ce (like sim- 
ply a ssuming the classical reference results bv iKolmogorovl 
Il94ll for energy spectra) is wrong or incomplete in most 
cases. Turbulence driving already emerged as a key issue 
in the context of compressible turbulence i n isothermal gas, 
relev a nt for the prob l em of star formation (jFederrath et al.l 
l20ld ). lKritsuk et al.l (|2010l ). on the other hand, argue that 
a universal scaling law should be observed at length scales 
that are sufficiently small compared to the forcing sc ale. This 
idea is supported by recent findings of Schmidt & F ederrathI 
(|201ll ). In this work, a subgrid scale model for supersonic 
turbulence indicates a scaling exponent for the unresolved 
turbulent pressure that is independent of the forcing. Nev- 
ertheless, the global statistics of turbulent pressure varies 
with the forcing. The physical conditions of the gas and the 
turbulence in the cosmological large-scale structure are obvi- 
ously different: here we are dealing mostly with subsonic or 
transonic flow. As this study shows, one needs the investiga- 
tion of the role of turbulence forcing and turbulent support 
in this regime, and their applications to cosmological simu- 
lations, in order to better understand the gas dynamics in 
the IGM. 
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APPENDIX A: DIFFERENT DEFINITIONS OF 
THE BARYON PHASES 

Throughout this work, the distinction between gas belong- 
ing to clusters or filaments has been based solely on baryon 
overdensity. We check here that the main results of this work 
do not depend on this definition by testing different charac- 
terisations of the baryons. In particular, we computed the 
evolution of the mass-weighted averages of the SGS turbu- 
lence energies of "clusters" and "fil aments" , a ccording to 
criteria based on temperature (e.g.. ICen & Ostriker 199^) 
and on total overdensity fe.g.. lVazza et al.l 2009a ). 

The temperature ranges defining the WHIM and IGM 



are lO'^ K < T < lO'^ K and T > lO'' K, respectively. In 
IVazza et al.l l|2009ah , the baryon phases are characterised by 
the total overdensity i5t = (p -I- /5dm)/pcr, where the suffix 
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Figure Al. Time evolution of the SGS turbulence energy for 
the phases labelled as 'clusters' and 'filaments' according to our 
definition based on baryon overdensity (solid lines) and on the 
criteria based on temperatu re (dashed lines) and on total over- 
density l lVazza et alj [2009al . dotted lines). The single lines are 
scaled by arbitrary factors, in order to ease the visualisation on 
the plot. 



'dm' indicates the dark matter density. In that work, fila- 
ments and sheets have overdensity 3 ^ St < 30, and the 
clusters are defined by St ^ 50. 

In Fig. lAll the evolution of et for the two baryon phases 
and the different definitions is shown. The definitions are 
not completely equivalent and therefore the average values 
differ from each other, in some cases up to an order of mag- 
nitude, but in Fig. lAll they are scaled for sake of clarity. 
The comparison shows that the trends in the time evolution 
of et presented and discussed in Section 13.21 are apparent 
also with other definitions of the baryon phases, and are not 
caused by the particular choice of using baryon overdensity 
phases. 
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